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SUMMARY 


This  report  describes  the  development  of  a  mathematical  model  for 
gunner's  tracking  performance  in  MTQ  Mode  II  tracking  task  which  is  a  linear 
time-varying  antiaircraft  artillery  system.  The  Luenberger  observer  theory 
is  used  to  design  the  gunner  model  which  is  composed  of  three  elements — a 
reduced-order  observer,  a  feedback  controller,  and  a  remnant  element.  An 
important  feature  of  the  model  is  that  its  structure  is  simple,  hence  the 
computer  simulation  of  man-in-the-loop  AAA  tracking  systems  using  the  gunner 
model  requires  only  a  short  execution  time.  A  parameter  identification  pro¬ 
gram  based  on  the  combined  least  squares  curve-fitting  method  and  the  modified 
Gauss  Newton  gradient  algorithm  is  developed  to  determine  parameters  of  the 
model  systematically.  Model  predictions  of  both  azimuth  and  elevation  tracking 
errors  for  several  target  flyby  and  maneuvering  trajectories  are  shown  to  be 
in  excellent  agreement  with  the  empirical  data  obtained  from  manned  AAA 
simulation  experiments  conducted  at  the  Aerospace  Medical  Research  Laboratory, 
Wright-Patterson  AFB,  Ohio.  It  is  concluded  that  the  antiaircraft  gunner 
model  based  on  the  observer  theory  can  be  accurately  and  efficiently  used 
to  study  AAA  weapon  effectiveness  and  aircraft  survivability. 
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PREFACE 


This  report  documents  a  study  performed  by  Systems  Research 
Laboratories,  Inc.  (SRL) ,  Dayton,  Ohio,  for  the  Aerospace  Medical  Research 
Laboratory  (AMRL) ,  Manned-Systems  Effectiveness  Division,  Manned  Threat 
Quantification  program.  This  work  was  performed  under  Contract 
F33615-76-C-5001.  The  Contract  Monitor  was  Mr.  Robert  E.  Van  Patten, 
and  the  Technical  Manager  was  Capt.  George  J.  Valentino.  The  SRL  Project 
Manager  was  Mr.  Charles  McKeag. 

The  authors  wish  to  extend  their  appreciation  to  Dr.  Carroll  N.  Day, 
Mr.  Walt  Summers  and  Dr.  Dan  Repperger  of  Manned-Systems  Effectiveness 
Division  of  the  Aerospace  Medical  Research  Laboratory,  WPAFB,  for  many 
valuable  discussions  and  helpful  comments. 
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SECTION  I 


INTRODUCTION 


The  Luenberger  observer  theory  [1]  has  been  applied  to  design  a  gunner 
model  by  the  authors  in  [2]  through  [8]  which  describes  the  gunner's 
tracking  response  in  antiaircraft  artillery  (AAA)  systems.  The  main  high¬ 
lights  of  the  model  are  a  simple  model  structure  and  accurate  predictions 
of  tracking  errors.  In  [4],  it  has  been  pointed  out  that  in  the  computer 
simulation  of  man-in-the-loop  AAA  tracking  systems  tremendous  computer 
execution  time  can  be  saved  by  using  the  human  operator  (or  gunner)  model 
based  on  the  observer  theory  rather  than  the  optimal  control  model  [9]. 

This  is  the  advantage  of  a  model  with  simple  structure.  The  designed  model 
has  been  applied  to  study  the  weapon  effectiveness  of  linear  time-invariant 
AAA  tracking  systems  at  the  Aerospace  Medical  Research  Laboratory  of  Wright - 
Patterson  AFB,  Ohio.  This  report  will  present  a  more  general  antiaircraft 
gunner  model  which  can  represent  the  gunner's  tracking  performance  in  the 
MTQ  System  which  is  a  linear  time-varying  AAA  tracking  system.  The  basic 
structure  of  the  model  contains  three  elements — a  reduced-order  observer,  a 
feedback  controller,  and  a  remnant  element.  However,  the  reduced-order 
observer  is  now  a  linear  time-varying  system.  So  the  antiaircraft  gunner 
model  is  also  a  linear  time-varying  system.  A  parameter  identification 
program  is  developed  to  determine  model  parameters  such  that  the  model 
output  can  represent  gunner's  tracking  function.  The  combined  least 
squares  curve-fitting  method  and  the  modified  Gauss  Newton  gradient 
iterative  algorithm  are  used  to  design  the  parameter  identification 
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program.  The  parameter  values  of  the  gunner  model  are  obtained  by 
iteratively  adjusting  their  values  until  the  model  predictions  of  azimuth 
and  elevation  tracking  errors  match  the  empirical  data  of  the  manned  MTQ 
simulation  experiments.  A  computer  simulation  program  for  the  closed-loop 
MTQ  tracking  task  is  also  developed.  Simulation  results  are  compared  with 
the  experimental  data.  Model  predictions  of  tracking  errors  are  shown  to 
be  accurate  representations  of  actual  gunner  tracking  data.  Many  figures 
showing  computer  simulation  results  for  various  flyby  and  maneuvering 
trajectories  are  also  included  in  this  report. 
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SECTION  II 


MTQ  TRACKING  SYSTEM 


Figure  1  shows  the  block  diagram  of  the  MTQ  tracking  system.  The 

visual  display  provides  the  information  of  tracking  error  e^,  on  a  two 

dimensional  visual  device  to  the  human  operator  (gunner) .  The  tracking 

error  e^  is  the  difference  between  the  target  angle  0  and  the  gunsigKt 

angle  Q  .  The  x-axis  and  the  y-axis  components  of  the  signal  on  the 

display  represent  the  azimuth  and  the  elevation  tracking  errors.  The 

function  of  the  human  operator  (gunner)  in  Figure  1  is  to  align  the  gunsight 

angle  G  to  the  target  angle  0  .  In  the  MTQ  system,  only  one  gunner  per- 
8 

forms  angle  tracking  task  (i.e.,  he  tracks  both  azimuth  and  elevation 
angles).  In  the  following  analysis,  the  MTQ  tracking  system  is  decomposed 
into  elevation  and  azimuth  tracking  systems  as  shown  in  Figures  2  and  3. 

The  gunner  models  in  these  figures  are  mathematical  models  representing 
gunner's  azimuth  and  elevation  tracking  responses  to  be  designed  in  the 
next  section.  In  Figure  3,  the  factor  cos(0  )  where  (0  )  is  the 

g  hL  g 

elevation  gunsight  output,  is  included  in  the  modeling  of  the  azimuth 
visual  display  for  some  optical  consideration  (see  Reference  [10]  for 
detail) .  Therefore,  the  elevation  and  the  azimuth  tracking  systems  are 
coupled  in  the  sense  that  the  elevation  gunsight  output  (0  )E^  is  used  as 
one  of  the  inputs  to  the  azimuth  tracking  system. 


The  input-output  relations  of  the  azimuth  and  the  elevation  gunsight 


systems  are  represented  by 


<V«  '  U,V 


(1) 


and 


<VeL 


1. 34u 


EL 


(2) 


Several  simulated  flyby  and  maneuvering  trajectories  of  the  target  aircraft 
described  in  [10]  are  used  as  inputs  to  the  AAA  tracking  system  for  simu¬ 
lation  purposes.  Now,  a  general  state  space  equation  of  the  azimuth  or  the 
elevation  AAA  tracking  system  including  the  gunsight  system  and  the  target 
motion  can  be  derived  in  the  following.  Let  us  introduce  state  variables 


Xil(t)  ^®T*i  ~  (<Vi  "  (0g^i 


x 
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A 

(t)  =  (6t)1 


where  i  =  AZ  or  EL  represents  azimuth  or  elevation  components  of  signals  in 

the  AAA  tracking  system.  For  example,  xA21  and  xA22  denote  the  azimuth 

components  of  tracking  error  and  target  angle  rate  respectively.  And  let 

T 

x i  -  [x  (t)  xi2(t)]  be  the  state  vector.  The  general  system  dynamic 

equation  and  the  measurement  equation  of  the  azimuth  or  elevation  AAA 
tracking  system  can  be  derived  by  using  Equations  (1)  and  (2)  and  Figures  2 
and  3. 
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x.  =  A.x,  +  B.U.  =  F. (Q  ) 
— i  i-Hc  i  i  1  T  i 


and 


yi  =  Ci^i 


(3) 

(4) 


where  i  =  AZ  or  EL, 


- 

r  - 

0  1 

0 

Fi  = 

0  0 

1 

-  - 

_  . 

with  bAzl  =  -1.28  and  bEL1  =  -1.34,  (5) 


A 

Ci  =  [cil  °]  With  CAZ1  =  C0S  (VEL  and  CEL1  =  (6) 


and  u^,  (0T)i  and  denote  azimuth  or  elevation  components  of  the 
gunner's  control  output,  the  target  acceleration  and  the  observed 
tracking  error  respectively.  The  above  matrices  A^  and  have  the 
same  values  for  azimuth  or  elevation  cases  but  the  vector  B.  in 

l 

Equation  (5)  shows  different  component  values.  Equation  (6)  shows  that 
CAZ1  is  time_varyinS  while  is  a  constant  equal  to  one.  Therefore, 

the  azimuth  MTQ  tracking  system  is  a  linear  time-varying  system  while 
the  elevation  MTQ  tracking  system  is  only  a  linear  time-invariant 
system.  In  the  next  section  gunner  models  will  be  designed  representing 
gunner's  azimuth  and  elevation  tracking  responses  corresponding  to 
equations  (3)  and  (4). 
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SECTION  III 


MTQ  TRACKER  MODEL  DEVELOPMENT 

The  authors  have  developed  an  antiaircraft  gunner  model  [2]  through 
[8]  based  on  the  Luenberger  reduced-order  observer  theory  [1].  The  basic 
structure  of  the  model  is  shown  in  Figure  A  consisting  of  three  main 
elements  -  a  reduced-order  observer,  a  state  variable  feedback  controller, 
and  a  remnant  element.  The  first  two  elements  contain  the  deterministic 
part  of  the  model  while  the  remnant  element  represents  the  random  (or 
stochastic)  part  of  the  model.  The  reduced-order  observer  gives  an 
estimate  of  those  unmeasurable  state  variables,  and  it  is  one  of  the 
characteristic  capabilities  of  a  human  tracking.  The  state  variable  feed¬ 
back  controller  represents  the  gunner's  control  function.  All  the  effects 
of  the  various  randomness  sources  due  to  human  psychophysical  limitations 
and  of  the  modeling  errors  are  equivalently  lumped  into  one  remnant  element. 
These  effects  include  the  observation  error,  the  neuromotor  noise,  target 
uncertainty  modeling  error,  tracking  error  due  to  tracking  difficulty  with 
respect  to  various  trajectories,  etc.  In  this  paper,  more  general  mathe¬ 
matical  equations  of  the  model  will  be  derived  corresponding  to  the  linear 
time-varying  tracking  system  described  by  Equations  (3)  and  (A). 

Since  the  gunner  doesn't  have  precise  information  about  the  target 
dynamics  (i.e. ,  the  so-called  human's  uncertainty  about  target  motion),  the 
term  representing  target  angle  acceleration,  0^,  in  Equation  (3)  will  not 
be  included  in  the  design  of  the  reduced-order  observer  of  the  gunner 


model.  Let  us  assume  that  the  gunner’s  understanding  (or  internal  model) 
about  the  tracking  system  is  described  by 


*i 


Ai-i  +  Biui 


and 


=  c<— 


i-i 


For  the  convenience  of  designing  a  reduced-order  observer,  let  us  introduce 


a  new  state  vectoi  xj 


x'  = 
-i 


cilxil 


12 


where  x^  and  x^  are  state  components  of  x^,  and  c^  is  the  first  element 
of  of  Equation  (6).  Then  the  equations  representing  the  gunner's 
internal  model  about  the  tracking  system  can  be  rewritten  as 


-i  i  i 


'i  “i 


(7) 


and 


L6 


y 


i 


CI*i 


where 


(8) 


- 

- 

_  “ 

cilcil_1 

cil 

cilbil 

Ai  = 

0 

0 

Bi  = 

0 

c;  =  [i  o] 

Now  Cj  is  a  constant  vector.  In  fact,  the  new  first  state  component  of  x| 
is  measurable.  Therefore,  only  the  second  state  component  x^  of  xj  needs 
to  be  estimated  in  order  to  implement  the  state  variable  feedback  controller. 
The  reduced-order  observer  theory  [1]  is  now  applied  to  give  an  estimate 
£^2  °f  xi2  using  Equations  (7)  and  (8) .  It  can  be  shown  that 
satisfies  the  following  differential  equation. 

*12  =  -kcil*i2  +  k*  -  k  ^il'1  y  -  k  Cilbiluc 


where  c^  and  b^  are  defined  in  Equations  (5)  and  (6),  k  is  the  observer 
gain,  and  u^  is  a  linear  feedback  control  law  of  the  form 


c  =  -[Y1  Y2] 


(9) 
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where  and  y2  are  two  feedback  control  gains.  The  values  of  k,  y^  and 
will  be  shown  later.  In  order  to  bypass  the  time  derivative  of  y  in 
Equation  (9),  let  us  introduce  a  new  variable  z. 


z  *  *12  -  ky(t) 


(10) 


then  the  observer  dynamics  can  be  expressed  by 


"kcilZ 


'  2 
k  c 


il 


kcilcil 


0- 


CilbilUc 


(ID 


Next,  the  control  output  u(t)  of  the  gunner  model  of  Figure  4  can  be 
expressed  by 


u ( t)  =  u  (t)  +  V ( t ) 
c 


+  v(t) 


(12) 


where  the  remnant  v(t)  is  modeled  as  a  white  noise  with  mean  zero  and  with 
covariance  function  as  follows 


E[v(t)  v ( t) ] 


a!  +  a2  °T  (t)  +  “3  ®T 


S(t-T) 


for  all  t  and  x 


(13) 


where  are  three  nonnegative  parameters  of  the  model  to  be 

£  •* 

determined,  0^,  and  0^,  are  estimated  target  angle  rate  and  acceleration 
respectively.  Note  that  by  definition  0T  is  xi2  and  can  be  obtained  from 
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Equation  (9)  or  Equation  (11),  §T  can  be  obtained  by  taking  a  first  order 
Taylor  series  expansion  of  8  .  Equations  (11)  through  (13)  are  the  mathe¬ 
matical  equations  of  a  general  gunner  model.  By  selecting  i  =  AZ  or  EL, 
this  model  represents  the  gunner's  response  in  azimuth  or  elevation  tracking 
task  respectively.  Next  the  gunner  model  equations  will  be  combined  with 
Equations  (3)  and  (4)  of  the  tracking  system  (which  includes  the  gunsight 
system  and  the  target  motion)  to  form  overall  equations  of  the  man-in-the- 
loop  AAA  tracking  system.  For  the  following  analysis,  the  state  components 
of  the  overall  system  are  selected  to  be 


X. 

— l 


(14) 


where  y^,  x^j  and  are  defined  in  Equations  (4),  (3)  and  (9)  and  k  is 
the  observer  gain.  Then  the  state  space  equation  of  the  overall  system  can 
be  described  by 

ii  ‘  Mi  +  £i  ( 5i)i  +  Vi  <15> 

where  A^,  F^,  and  matrices  are 


CilCil  +  Cilk  "  bilCil(Yl  +  Y2k)  Cil  _bilCilY2  bilCilY2 
2  1 

"Cilk  -kdilCU  +  bilkcil  (Y1  +  Y2k)  -cilk  +  bilY2kcil  ~bilY2kcil 
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"o' 

rbiicn  1 

1 

.  D.  = 

u 

’  -i 

il  il 

_1_ 

. 'bilkcil. 

and  v  is  the  remnant  element  in  which  the  parameters  a^,  »  and  of  its 

covariance  function  may  have  different  values  for  azimuth  and  elevation 
tracking  cases.  The  parameters  of  the  gunner  model,  i.e.,  k,  y^,  Y2*  a^»  a2> 
and  of  Equations  (11),  (12),  and  (13)  will  be  determined  in  the  next 
section. 
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Section  IV 


COMBINED  CURVE-FITTING  IDENTIFICATION  PROGRAM 


Before  the  gunner  model  based  in  the  observer  theory  can  be  used  to  describe 
gunner's  tracking  response,  the  parameters  k,  y^,  ct^ ,  c^,  ,  of 

Equations  (9)  through  (13)  should  be  determined  first.  The  authors  have 
developed  a  curve-fitting  parameter  identification  program  using  the  least 
squares  method  [11]  and  the  Gauss  Newton  gradient  method  [12],  Model 
predictions  of  the  ensemble  mean  and  the  ensemble  standard  deviation  of 
tracking  errors  will  be  used  in  the  curve-fitting  program.  The  MTQ 
elevation  tracking  system  will  be  used  an  as  example.  The  corresponding 
equations  are  derived  in  terms  of  model  parameters  in  the  following. 


Taking  expectation  values  of  both  sides  of  Equation  (15),  we  have 


EL 


*ELXEL 


-el(°t)el 


(16) 


where  X_T  =  E  [X  ] .  The  first  component  x  of  X  is  the  model  prediction 
EL  EL  1  IjI-* 

of  ensemble  mean  of  tracking  error.  By  solving  vector  differential 
Equation  (16),  x^  can  be  expressed  in  terms  of  model  parameters.  (x^(0) 
is  assumed  to  be  zero.) 


x1(t)  = 


k  +  1,34t1  +  1-34k>2  1.34y  (t 

Jo  [l.34Yl  (k  +  1.34Yl)  6 


t) 


+ 
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1. 34y2 
k  +  1. 34y^ 


-k(t  -  x) 
e 


1  +  1.34y2 

1.34Yi 


6t(t)  dx 


(17) 


Furthermore,  let  P(t)  =  E  jjXgj/t)  -  ^(t))  (x^Ct)  -  Xg^Ct))  TJ  ,  and  it  can 
be  shown  that  the  first  diagonal  element  of  P  is  the  square  of  the  ensemble 

standard  deviation  of  tracking  errors.  It  can  be  shown  that  P(t)  satisfies 
the  following  covariance  equation. 


P  -  AelP  -  P  ^EL  -  5EL  (  «1  +  “2  VC>  +  “3  V'>)  2*  EL  <18) 

By  solving  this  matrix  differential  equation  and  with  the  assumption  p  (0)  =  0 
the  first  diagonal  element  p^  (t)  of  P(t)  can  be  obtained. 


f  1.342  • 

\  -1 

O 

I  k  +  1.34y1 

1.34kY2  e 


■k(t  -  x) 


\  !•  3^Y-|  (t  -  t) 

(k  +  1.34y  +  1.34kY2)  e 


(o^  +  a2  0T(x)  +  a3  eT(x)dx) 


(19) 


Let 


s(t) 


Pll(t) 


1/2 


(20) 


22 


then  s(t)  denotes  the  model  prediction  of  ensemble  standard  deviation  of 
tracking  errors.  Equations  (17)  and  (20)  are  ensemble  mean  and  ensemble 
standard  deviation  of  tracking  errors  which  are  explicit  functions  of  model 
parameters. 

The  parameters  of  the  antiaircraft  gunner  model  will  be  determined  system¬ 
atically  and  simultaneously  by  a  combined  least  squares  curve-fitting 
identification  program.  The  reference  curves  to  be  fitted  in  the  curve¬ 
fitting  program  are  obtained  from  empirical  data  of  manned  AAA  simulation 
experiments  conducted  at  the  Aerospace  Medical  Research  Laboratory  of 
Wright -Patterson  AFB,  Ohio.  Several  simulated  flyby  and  maneuvering  air¬ 
craft  trajectories  of  35  seconds  duration  are  used  as  target  trajectories 
for  the  above  experiments.  Let  x|(t)  and  s'(t)  be  the  reference  empirical 
mean  and  standard  deviation  of  tracking  errors  which  are  obtained  by 
averaging  the  results  of  twenty-five  experimental  simulation  runs  with  the 
same  target  trajectory  and  the  same  subject.  Now  the  criterion  function  J 
of  the  combined  least  squares  curve-fitting  program  is  defined  as 


where  t^  is  the  tracking  duration  (equal  to  35  seconds  in  this  case).  x| 
and  s'  are  empirical  reference  data,  is  the  parameter  vector  of  the 
gunner  model. 


(21) 
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k  Y1  y2  al  a2  a3 


x^  and  s  are  functions  of  time  and  parameter  vector  a^  as  shown  in  Equations 
(17)  through  (20).  c  is  a  positive  weighting  factor  selected  to  be  one  in 
this  case.  The  criteron  function  J  of  Equation  (21)  can  be  considered  as  a 
combination  of  two  integrations.  One  of  them  is  an  integration  of  the 
square  of  the  error  between  the  empirical  mean  tracking  error  x|  and  the 
model  prediction  of  the  ensemble  mean  tracking  error  x^.  The  other  is  a 
similar  integration  related  to  the  standard  deviation  of  tracking  errors. 
Therefore,  minimizing  J  with  respect  to  the  parameter  vector  a_  is  called  a 
combined  least  squares  curve-fitting  method.  We  would  like  to  point  out 
that  in  this  method  model  parameters  will  be  determined  such  that  empirical 
mean  and  standard  deviation  of  tracking  errors  are  fitted  simultaneously  by 
the  corresponding  model  prediction  functions  x^  and  s.  Now  the  parameter 
identification  task  becomes  a  minimization  problem,  i.e.,  to  find  values  of 
the  parameter  vector  ji  which  minimize  the  criterion  function  J(a)  of 
Equation  (21).  The  modified  Gauss  Newton  iterative  gradient  method  is 
derived  in  Appendix  A  which  iteratively  adjusts  parameter  values  to  mini¬ 
mize  the  criterion  function  J.  This  iterative  process  will  continue  until 
the  increments  are  smaller  than  a  preassigned  lower  bound.  A  computer 
curve-fitting  program  is  developed  using  the  above  mentioned  methods  and 
procedures.  The  highlights  of  the  identification  procedure  are  a  fast 
method,  accurate  values  and  a  systematic  approach.  Obviously,  the  parameter 
identification  done  by  a  computer  is  much  faster  than  manual  tuning.  The 
combined  least  squares  curve-fitting  minimization  procedure  can  definitely 
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provide  more  accurate  parameter  values  than  the  trial-and-error  approach. 
Finally,  it  is  a  systematic  procedure  to  determine  parameter  values  and  can 
be  applied  to  gunner  models  describing  control  responses  for  various  AAA 
tracking  systems.  The  parameter  values  of  the  gunner  model  determined  by 
this  identification  procedure  for  both  elevation  and  azimuth  tracking  tasks 
are  listed  in  the  following  table. 


Parameter 


Observer 

Gain 


Gunner 
Model  For 


Elevation 

Tracking 


Azimuth 


Controller 

Gains 

Coefficients  of 

Remnant  Covariance 

Function 

y  1 

Y2 

«1 

a2 

a  3 

-1.99 

-.745 

.0000094 

.025 

.068 

-3.51 

-.762 

.0000363 

- 1 

.00614 

.0117 

We  would  like  to  point  out  that  the  above  parameter  values  are  obtained 
consistently  from  the  combined  curve-fitting  parameter  identification 
program  with  various  initial  guesses.  In  addition,  the  parameter  values  of 
the  gunner  model  depend  on  the  dynamics  of  AAA  gunsight  systems.  For 
different  AAA  gunsight  dynamic  systems,  the  combined  curve-fitting  program 
will  determine  different  parameter  values  for  the  gunner  model. 


Section  V 


COMPUTER  SIMULATION  RESULTS 


Once  the  parameters  of  the  gunner  model  are  determined,  the  gunner  model 
is  ready  to  be  used  in  computer  simulation  of  the  AAA  closed  loop  tracking 
system  to  describe  gunner's  tracking  performance.  These  parameter  values 
can  be  substituted  into  the  elements  of  matrices  A £,  F^,  and  of 
Equation  (15)  which  is  a  mathematical  dynamic  model  of  the  overall 
closed-loop  AAA  tracking  system.  In  order  to  find  a  solution  of 
Equation  (15)  without  using  convolution  integrations,  it  can  be 
discretized  to  be 


where 


x 

— n+1 


x  +  r 


0T  „  + 
T,n 


(22) 


with  tn+1  =  (n+1)  At  and  At  =  .066  seconds. 


(23) 


r 


2 


0  =  6  (t  )  and 

T,n  T  n 


vr  is  a  random  sequence  with  the  following  properties 
E  [v„]  -  0, 

E  [(’«)  (vn)T]  '  T£  (“l  +  K  K)  +  “3  ('„)) 


(24) 
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Taking  expectation  values  of  both  sides  of  Equation  (22),  we  get 


n+1 


X  + 
n 


F1  9T,n 


(25) 


Let  P  denote  the  covariance  of  X  then  it  can  be  shown  that  P 
n+1  —n+1 

satisfies  the  following  matrix  difference  equation. 


n+1 


n+1 


-  ♦  p  ♦  +  r„ 


JL 

At 


(°1 


+  a  0_ 


+  «3  eT 


(*.)) 


(26) 


Then  the  first  element  of  X  of  Equation  (25)  and  the  square  root  of 
the  first  diagonal  element  of  the  matrix  Pn+1  of  Equation  (26)  are  the 
model  predictions  of  the  ensemble  mean  and  the  ensemble  standard  deviation 
of  tracking  errors  respectively.  Computer  programs  simulating  MTQ  closed 
loop  tracking  system  are  developed  using  the  above  recursive  equations 
(25)  and  (26).  The  inputs  to  these  programs  are  motion  trajectories  of 
target  aircraft.  Six  flyby  and  maneuvering  trajectories  as  shown  in 
Figure  5  are  used  in  this  study.  The  antiaircraft  weapon  is  located  at 
the  origin  of  local  horizontal  x-y  plane  and  the  z  axis  represents  the 
altitude.  The  increment  of  each  of  the  three  axes  is  1000  ft.  A  detailed 
description  of  the  characteristics  of  these  trajectories  can  be  found  in 
[10] .  The  output  of  the  computer  simulation  program  are  ensemble  mean 
and  ensemble  standard  deviation  of  tracking  errors  which  can  be  used  to 
predict  the  gunner's  tracking  characteristics.  Computer  simulation 
results  of  the  MTQ  tracking  system  show  that  model  predictions  match 
well  with  the  empirical  data  for  both  mean  and  standard  deviation  of 
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FIGURE  5. FLYBY  AND  MANEUVERING  TARGET  TRAJECTORIES 


tracking  errors.  Simulation  results  are  shown  in  Figures  6  through  29  for 
the  trajectories  of  Figure  5.  There  are  two  curves  on  each  of  these 
figures.  The  solid  curve  is  the  empirical  data  which  is  obtained  by 
averaging  the  results  of  twenty-five  experimental  simulation  runs  of 
35  seconds  duration.  The  dashed  curve  shows  the  ensemble  value  of  the 
model  prediction  of  the  tracking  error.  Figures  6  through  9  show  the 
results  of  elevation  mean,  elevation  standard  deviation,  azimuth  mean, 
and  azimuth  standard  deviation  of  tracking  errors  for  the  maneuvering 
trajectory  T1  of  Figure  5.  All  of  these  four  figures  show  that  the  MTQ 
tracker  model  can  provide  accuracy  prediction  of  the  corresponding 
empirical  data.  Next,  Figures  10  through  13  show  the  similar  results  for 
the  flyby  trajectory  T2  of  Figure  5.  Again  three  figures  indicate  the 
designed  MTQ  tracker  model  generating  model  predictions  in  excellent 
agreement  with  the  empirical  data  for  a  flyby  trajectory  too.  We  would 
like  to  point  out  that  the  same  set  of  parameter  values  determined  by 
the  combined  curve-fitting  identification  program  in  the  last  section  can 
be  used  to  generate  accurate  model  predictions  for  all  six  flyby  and 
maneuvering  trajectories  of  Figure  5.  Therefore,  it  is  concluded  that 
for  a  given  gunsight  weapon  system  (e.g.,  MTQ)  a  same  set  of  parameter 
values  can  be  used  in  the  observer  model  to  predict  human  operator's 
tracking  errors  for  all  realistic  simulated  target  trajectories.  Hence, 
the  observer  model  (i.e.,  MTQ  tracker  model)  is  a  predictive  model. 
Similarly,  Figures  14-17,  18-21,  22  -  25,  and  26  -  29  show  the 
corresponding  simulation  results  for  trajectories  T3,  T4,  T5,  and  T6. 

All  these  results  verify  that  the  designed  observer  model  can  describe 


human's  tracking  response  accurately.  The  names  of  these  six  trajectories 
in  [10]  are  listed  as  follows  for  reference. 

T1  Zigzag 

T2  2x2  Flyby 

T3  Fairpass 

T4  RECON 

T5  Weapon  Delivery 

T6  4  x  .3  Flyby 


A  comparison  between  this  model  prediction  and  the  optimal  control  model 
prediction  [9]  is  shown  in  Figure  30.  The  dotted  curve  in  Figure  30  is 
the  optimal  control  model  prediction  of  the  tracking  error.  Obviously, 
it  shows  that  the  antiaircraft  gunner  model  can  represent  the  gunner's 
response  as  well  as  that  by  the  optimal  control  model.  However,  the 
computer  execution  time  by  simulating  the  linear  time-varying  MTQ  gun 
system  using  the  gunner  model  is  less  than  15%  of  that  using  the  optimal 
control  model. 
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Section  VI 


CONCLUSION 


This  report  presents  some  extended  work  on  the  development  of  an  anti¬ 
aircraft  gunner  model  based  on  the  observer  theory  for  the  MTQ  system. 

In  [2]  through  [8],  the  authors  developed  a  gunner  model  for  the  study  of 
weapon  effectiveness  of  linear  time- invariant  AAA  systems.  The  more 
general  antiaircraft  gunner  model  developed  in  this  report  can  be  as 
well  applied  to  linear  time-varying  AAA  systems  (e.g.,  MTQ  system). 

The  structure  of  the  gunner  model  contains  three  elements  -  a  reduced-order 
ovserver,  a  feedback  controller,  and  a  remnant  element.  Here  the  reduced- 
order  observer  is  designed  as  a  linear  time-varying  system.  A  computer 
simulation  program  is  developed  with  the  designed  model  describing  the 
gunner  tracking  performance  for  a  linear  time-varying  AAA  tracking  system. 
Computer  simulation  results  are  in  excellent  agreement  with  the  manned 
AAA  simulation  empirical  data  for  several  flyby  and  maneuvering  trajec¬ 
tories.  It  is  also  shown  that  the  gunner  model  can  predict  tracking 
error  as  accurately  as  the  optimal  control  model.  However,  the  computer 
execution  time  of  the  MTQ  closed  loop  tracking  system  simulation  utilizing 
the  gunner  model  is  much  shorter  than  that  using  the  optimal  control 
model.  All  these  results  verify  that  the  designed  gunner  model  is  an 
accurate  and  efficient  model  describing  the  gunner's  compensatory 
tracking  characteristics  for  various  linear  time-varying  or  time 
Invariant  AAA  weapon  systems. 


The  six  model  parameters — an  observer  gain,  two  controller  gains,  and 
three  coefficients  of  the  remnant  covariance  function — are  determined 
by  the  combined  least  squares  curve-fitting  identification  program. 

This  identification  program  provides  an  efficient  model  validation  method. 
The  designed  gunner  model  and  the  parameter  identification  procedure  have 
been  successfully  applied  to  study  the  problems  of  the  aircraft  surviv¬ 
ability  and  air  defense  weapon  effectiveness  at  the  Aerospace  Medical 
Research  Laboratory,  Wright-Patterson  AFB,  Ohio.  All  the  results  of  the 
computer  simulation  of  closed-loop  AAA  tracking  systems  show  that  model 
parameters  can  be  determined  accurately  by  the  curve-fitting  program  and 
the  gunner  model  is  a  reliable  representation  of  actual  gunner's 
tracking  performance. 
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APPENDIX  A 


Derivation  of  the  Iterative  Algorithm  of  Modified  Gauss  Newton  Method: 
The  criterion  function  J(a)  of  Equation  (21)  is  rewritten  here, 


J(a) 


[(*;«> 


-  xj(t,a))2 


+  c  (  s' (t) 


dt 


Taking  first  order  Taylor  series  expansions  of  Xj  and  s  with  respect  to 
a  certain  initial  guess  a^,  we  have 


+ 


— 

2 

/  (t,an) 

c 

»'“>  -  •  (t-lo)  -  • 

- 

_ 

4 

Next,  the  partial  derivative  of  J  with  respect  to  a  can  be  found  as 
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2(s'(t)  -  s  (t ,_ao)^ 


3s  t*aQ 


3s(t,ao)  3s  (c,ao) 

2  -  (a  'A>)  - dt 


Assume  that  a*  is  a  minimal,  then 


3 J(a*) 


After  discretization,  Equations  (Al)  and  (A2)  can  be  rearranged  to  be 


£  2 


•srft).-(sni 


j  ”2  (xl  (t)  ~  xl)  \  3a 


9?.  \T 


-2  (!■  -  =)  (|f 


Equation  (A3)  can  be  extended  to  a  more  general  form, 


i  =  0,  1,  2,  3, 


This  is  the  modified  Gauss  Newton  iterative  gradiant  algorithm. 


« 
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APPENDIX  B 

PRGRAM  LISTING  OF  PARAMETER  IDENTIFICATION  PROCEDURES-AZIMUTH  CASE 


H8K,T6»C  r*Ol(IO.CMtlO»*0  ,6TCS8.  k7  **  3  3  s  ,  *0 ~  r8»<- gfrfrfc - 

FTN. 

COMMENT,  HG  peck—**** - 

ATTACH, TAP£l,CAYiC,GY=2. 

AT-TACH,  IMCL  ,  TC=t  ABt' A^-Y-,  SN=A4>0 ,  Mis=l « - 

LIBRARY, I KSL. 

LUC, - 

PROGRAM  AZFIT  (INPUT .OUTPUT  ,'f  APE1) 

- C&HKCN/ VA../0T  ,N1,  L2  ,  IP,  II  Ml-,4  I  ,  kF - 

COMMON’/ fAri/F 2 <  El 2, 5)  ,»>  <6  ) ,  -.A  <6)  ,£  (  512, 14) 

- C0H*CN/»m^F44-20  48  >,^<C)  ,P*iJ0  <  &Ar24 - 

PATA  1P,WF/6,1../ 

C— US  £  -EV-SrfY-  V 1 D’  0  A  >  A — B  c;uT - 

c  EAD* ,11 

- RlaMAl - 

N2=1+1C A9/I1 

- IF-4U2,  Cl  .-E-l  2  >  *’4-2*  512 - 

0T=.  <533*11 

- 3KC1-  ,  Al-,3>.,r.Y--fUT  »  U.--r3 w U 2 - ,-t»2 - 

DO  15  1*1,2 

- ..ici^a - 

no  10  <1=1,1050 

- £6-0=:  Zfr - - - 

P£AS<1,*>  Z1,72,Z3,Z4,?5,7E 

- JCttaJO*!  1-1 - 

IF(MCG(«,I1>.ME,u)  GO  “.0  1C 

- X^laX^lAl - - - - 

IF  (IC1.GY.512)  GO  TO  10 

- IF«,e0.2>  GC  £G--§> - 

F2  (IC1»1)=Z1 

- r?U  C1,2»=Z2 - 

F  2 <ICl, 2) =Z6 

- B  (I-Cl  ,1  >*Z3 - - 

PHI'3(lCl)=(Z6-Z6u  )/,0  33 

- CO-TC  13 - - - - 

5  72  =  Z2/(U,CO.*OOS<BUC1,1))) 

- 7faaaw<i,ua.*r.a-HF. cum  >  > > - 

F2(ICI,4)=Z2 

_ F2UC1,  EU74 _ 

10  CONI INUE 

- IS  CONTINUE _ 

C  INPUT  "i”  10  GC  Gi-.  "0"  TO  SrCF 

- 20  REAUf  »IFf> _ _ _ 

IF  (IFG, t0,0 )  LT0F 

C  MA<F  INITIAL  ..GUCSS  AND  iU£ Ul  ID'.  A^iCN  .LIMITS - 

C  INPUT  LlMl,  L1M2 ,  V  EE 

- N6AO»,llKltLlH2-r€fc - 

PPINT*,5hLIM1  =  ,UM1,5HL:n2=,1  IM?,  2HEE=,  EF 

o  INPUT  1ST  r.tir JU-K&AXau  4iPt.-«iT 

READ*,*.  A 

- PA  -IN  I  ♦  ,  4  HI  ST  C.!X£^  .  ,  A  A _ _ _ 

C  COMPUTE  INITIAL  GC “T  FUNCTION 


PRINT*, 6HJHIN1*,RJ2 

- CALL  LOOP  <R«12) - 

GO  TO  20 

- EWJ - 

SUBROUTINE  LCOPCRJ2) 

C  IT ERA TICK  PROCESS - 

DIMENSION  R<6) ,0(6,6) ,QI (6,6) , HCkK ( 6 ) ,DDL (6 ) 

COMMON/ VAR/0T,Ni,  N2, IP, LIM1,LIM2, EE,  WF 
COMMON/ HAT  i/F  2  (  ?12,5>  ,A  (6)  ,  A  A  (6)  ,Et  512,10 

- GO  HH<m/  MAT  2  /f  3  *20  40-WP  <  fc )  ,  PMl-B4Gt 2 ) - 

NCT=MCT=0  $  D£l=l.P 

- 5-  00-19  1*1,  IP - 

PCII=0. 

- D«I>*P, - 

DO  10  J=1,IP 

- IB  Q II ,  J)  *3  « - 

00  20  rf=Nl,N2 

- FOIF»-Gi » <  F  2 i  K  Hr>  -6  f  K ,  7 ) ) - 

SDIF=-2.**F<MF2 (K,5)-B(K,14)  ) 

- BO  *1  *  IP - 

R(I)=R<IMB(K,I  l*EDIF*B  <<»I  f7  )*3CIF 

- DO  15  J«1,IP - 

15  Q(X,J)*Q(1, J)*2«*B(K,1)*B  (K,  J)*2«*WF*B(K,i4>7)*B(K,  J*7) 


C  USE  IMSL  LIE  ROUTINE  "LINVIF*'  TO  COMPUTE  TtiVll  St  OF  G 

aiii  i  «  k  iij  ^  r  r  r  r*  a >  ha'  u 


XF(1ER.EQ.129)  PrtINf* ,4HIEk=» I: 


30  DO  32  1=1, IF 


RETURN 


40  00  45  1*1, IP 


45  IF (000 (I ) «GT. £E )  GO  TO  5* 


50  MCT=MCT4i 


PRINT**  7HERRCR  2,lilH  NGT  =  ,NCT »  5H  MCT=,HCT 


PRINT*  »  5H  RJ1=,RJ1,5H  RJ2=,KJ2 


55  PJ2=PJl 


65  AA  (I)  =MI ) 


65  PRINT  70,A,NCT,XCT 


PRINT*  »  JHAA* »  aA 
PRINT*, 3H  D=,C 
PRINT*. 5NJMIN=, RJ1 


SUBROUTINE  CCEff(N,<J> 


01  MENTION  THC(512>  ,  fHD0(512>  ,THl  k  151 2 )  ,  I  ECO  K  <  5 12  ) 


C0KMCN/MATi/F2( El  2, 5) ,A (6), Am (6) ,5 ( 512,14) 


REAL  K 


ALFl  =  W(t|)  S  ALF2*M5>  $  AlF3*M6> 


T= (1-1) *07 


Cl=CCS(F2(I,3)) 


TwnTimZI 


10  F3 ( I) 


CO  2&  2=1, N2 


IFCI.EQ.l)  GC  TO  20 


20  CONTINUE 


CC  JO  2=1, N2 


9wm=E2(2,2) 


t\ N  (F2(I,3>  )*FMir  <I> 

62 

IF  (ARG1.LT.  2fcu.)  Sl  =  EXP  (-AkC-i  ) 


If  irlfrJlm 


X3»-T*X2-1. 


miaiu  ifrieki 


u  I  ■At  •  4  Fl  C  m.UdUft  Liffll  M IX  MLFUmfXrl 


CALL  VCCNV0<F3,BW,N2,N2,IHK> 


T=(I-i)*CT 


B  (I *1 ) =F 3 (I ) *Dt 


IKrTHik'flJ'XIVlLlll 


C2=-TAN  IF2 (I #  3) )*FHID«I) 


k  Fltri: i  WiAJK  Lirt  l.bl 


X2*X1*K*C1 


i  in  miUtf  Uf  n  n  ifiiiiii  urf  Fit  I 


X4*  (!•  28*CD  #*2*G£M2/X24*2 


IrC  MnLicJ 


ARG2*X1*T 


IF<ARG1.LT.200.  »  Si=EXP(-ARGi) 


HJfl.ftrJltt  fOMiMM.tiawatraiWI 


40  F3(I»  =  X3*  <S2*<Xl*T-l.m.)  +  X4* <S2*  <-X2*T +1 . )-Si> 


TiMirma 


00  50  I=i, N2 


t  rv4rllkiil  i 


BH(I)=F2(I, 2) 


C1=CCS  (F2 (1 1 3) ) 


riixjiWLfiiJaiuimi 


Xi*C2«-i.2a*Ci*GAMi 

X2=X1*<»C1 

X3si.2«*Cl/Xi 


ARGl*Cl*K*T 


Si*0.  $  S2«0. 


IF<ARG2.GT,-2C0«>  S2=EXP(ARG2 


CALL  VC0NV0(F3,8W,N2,N2,IWK) 


B(I.3)*F3III*OT 


9<I,5>=0. 


T*<X-1>*0T 


Cl*C0S(F2(X» 3)) 
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X1=C2*1.28*G1*GAM1 

— n-e..)(.i-+*«ci - 

X3=CC1+1.28*C1*GAM2)/X1 

- XV»  -*-r2^*C  1  ♦  0*  M  2/  »-2 - 

ARGi=Ci*X*T 

- Af?03»*-1«T - 

Sl=a.  $  52=0 • 

- I-f  (APGl«i.T,ge-(  .  )  -S-1  =  EXP(-A»  01  ) - 

IF  (ARG2.GT.-2tO.)  S2=fcXPC/RG2) 

68  F3  fI4« X  3* <a  2*1 1 > ♦ >4* (£  g-Sl ) - 

CALL  VCC.WxO<F3t3W,N2,N2,*WK) 

- DO  70  I«1>N2 - 

6  CI,7>  =F3(I)*LT 

- F»<i-1H»0T - - 

8W ( I) =F  2 (1 , 2 ) 

- 01«0e£(-Fg(li3-)) - 

AF.Gl=ci*x*r 

- Si-0. - 

IFCARGi.LT. 2«C)  Sl=cXP(-A^Gl) 

70  f 8<1)=£1«T»S1 - 

CALL  VCCNVOCFSt^W.NatMg.IWO 

- DO  80  I-1.N2 - 

THDK(I>=F3(I)*DT 

- If  <1.60.1 )  GC  TO . AS - 

THDDKCJ  )  =  (THC<(I)  -VHCK(I-l)  )/:■'( 

80  OOHT-IXOg - 

THOOKC1)=7HOOK(2) 

- DO  89  -1«1,N2 - - 

T=  <  1-1 )  *CT 

- ORa)»ALFl»ALF3*T-HP<I)«»g-+AlF3»TFHOD  <l)»*g - 

85  F3CI)=FT(F2(If3)fFHir(t) ,Y, GAH1.GAM2. K) 

- CALL  VCONVOCF-3,  BW , *g ,Kg , I KK ) - 

no  90  I =1 f N2 

- ^<1 =SGRT  <^3  <14*0T  ) - 

T*<I-1)*DT 

- W(l)  . ■  ALF 3*T H>4* > HW4t )  +  2  . » A  l  F  g «T HOl-4 i  )  •')  Hi- L  A  i-f 4 

90  F3CI)*FT(F2CI,3),FHID(I)  ,T  ,  C  GAH2  ,  K) 

- CALL  VCOHVOCF3.6W-rN3»Ng,  I-WK4 - 

DO  100  1=1 » N2 

— T-a-i>«ci - 

8<It8)=F3<I)*DT 

- BWCl>-ALFl»ALFg»T»*r.<l)»<r2*A4iF3.»rHOOU.lt.f.2 - 

Cl=COS  « F2 (X»  3) ) 

- C2  =  -TAMCF2<I,3)  )*PHI0(1) - 

Xl=C2*1.28*Ci»GAMl 

- X2«xl»i<*Cl - 

S1  =  0.  $  S2=0  • 

ARGi*Cl«K*T 

ARG2=xm 

— xr carci.lt. oco. >  si=exPc»AP.6i) - 

IFCARG2.G1.-2ta.)  S2  =  EXP  (AF.G2 ) 

- XI«1.2MCl»K«CA*aV-*2 - 

XA»tl.28*Cl) **?*2 .*<<!• ♦X3)*S2-X3*S1) 


64 


X2=Xi+<*Cl 


ARGi*Cl*K*T 


IF<ARGi.LT.2Pt .  )  S1»EXP (-ARG1) 


X3sl.26*Ci*K*GAM2/X2 


X5*2.*Ul.  +  X3)*s2-X3*Si)*(i.26*ci)*»3 


CALL  VCCN\0(F3| EW»N2»N2»IWK) 


T=a-i)*CT 


B(I,9)=F3(I)*0T 


C2*-TAN(F2(I,3) ) *  FHID  (I) 


X2sXl+K*Cl 


ARG1=C1«K*T 


IF  (ARGi«LT«2(ft*  )  $1*EXP  (-AKG1 ) 


X3=i.2fl*Cl*K/X2 


X5*<1.28*C1>**Z*Z.M<1.*X4)*$2-X4*S1> 


CALL  VCCNVC<F3tBW,N2,N2,IWK> 


T=(I-1)*DT 


B(I,1P)=F3(I)*QT 


CALL  VGCNV0(F3,PH,N2,N2tIWK) 


T*  C I- 1 1 *CT 


CALL  VCCNV0(F3,?W,S2,»:2,IKK) 


B(I»i2)*F3<I)*DT 

150  F3-(I)-=FMF24-I»3)  »PHIB<I)  »T » GAM1»£AM2 »K) - 

CALL  VCCN V0( F 3»  BK  » K2  |N'2 *  I WK) 

- DO  16Q  - 

160  B  < I »13) =F3 ( I > *0T 
DO  165  J=6, 13 
00  165  1*1 »  N2 

163-fr4*r>H  »f  5*B-C  I  t-sH/EUiIA) - 

RJ*0. 

- DO  17ft  - - 

170  RJ»RJ* (F2(I»4)*B<I*7) }**2  *WF* (F2(If5)-B(l»14))**2 

- B+3 - 

FUNCTION  FT(ThT ,PH0,T.GA«1,G4M2,K> 

- ei-cos  <w> - 

C2=-TAN (ThT)*PH0 

- XI  =c?+  l-.-g  a  ♦€  1  »6  A  M  l - 

X2=X1+K*C1 

- s,.  -qfr-ft, - 

ARG1*C1*K*T 

- A-RG3*X1*? - 

TF<AKG1.LT. 2Rt.  )  S1  =  EXP(-AFG1) 

- If  f  A&GM1  w-  BP- Dr  >  - -S  2°  ttt?  <  Afr-6  5  i - 


X3=1.28*C1*GAM2*K/X2 


APPENDIX  C 

PROGRAM  LISTING  OF  PARAMETER  IDENTIFICATION  PROCEDURES-ELEVATION  CASE 


H£LP,CY  =  3,IC  =  L7AG23ti 

W8F,T800, 10400,  CM700D0 , STCSP .  L74C2  38,GLra$,^tE»imir -  - 

FTN,L  =  0 • 

ATTACH.IMSL, io=library,$n=aso.  - - - 

LIBRARY, I MSL. 

ATTACH, T,DAT10,CY=1,WR=1.  - - - 

C3PY,T,TAPtl. 

LGO.  - 

PROGRAM  FIF< INPUT, OUTf UT.TAPtl ) 

C0MM0N/VAK./LIM1,LIM2,FE,IP»H,N,D£L,'N1  *C£ - - - 

C0MM0N/MAT/F2<512,h) , A ( 6 ) , A A ( 6 ) , B (51 2 ,9 > , 0 ( o > , F 3 ( 51 2 , 4) 

REWIND  1  - - 

C  CURVE  FIT  PROGRAM  FOR  ROOMS  -  MTQ  -  ELEVATION 

C  parameters  identified  are*  - - - 

C  A  ( 1)  =<GA1N  A(2)=GAMi-Al  A<3)=GAMMA2 

C  A  (4)  =A  LPHA1  f  CL;)  =  ALPHA  2  A  <  6 )  a  ALPHAS  ‘  ' -  . 

C  SAMPLES  OF  DATA  ARE  S. LOOTED  AT  EVERY  II  POINTS 

C  TIME  STEP  =  II*. 0  33  -  - - ~  "  - . 

C  P IKST  THREE  SECONCS  OF  DATA  APE  CONSIDERED  AOUISITION  AND  NOT  USE 

C  CE=  STANDARD  DEVIATION  WEIGHTING  FACTOR  TTTE*^I7T  - ' 

PRINT*, 30H  PIC<  tVERY  --  SAMPLES 

REA  J*,  II  -  — . -  -  . . - 

N=1+1U49/I1 

IF (N, GT«  5 12 ) N=512  . -  - -  ' 

H  =  . 033*11 

Nl  =  90/Il  - ~~  - - 

IP=fc. 

PRINT*  ,3Hil=  ,11, 4HOFL=,H,3H  N  =  ,  N ,  3H  N 1  = ,  N 1 ,  &Ti  T  P  = ,  I P  - 

C  READ  IN  TRAJECTORY  INFORMATION  AND  EMPIRICAL  LATA 

C  ALL  DATA  IN  rAOIANS 

HO  6  1=1,2 

DO  -j  <1=1,10  50 

2  «EAD(1 ,*) Zl, Z2,73,Zn,/5,Zb 

<<=<:+ii-i 

IF  (MOD  C<K,I1>«NL*0)GO  TO  6  . ‘  " 

4  1C1=IC1*1 

IF  ( ID  1  ,GT  ,512)  SO  TO  6  —  - - - - - -  - - 

IF  (I, E  0,1)G0  TO  £> 

ZU=Z3/1000.  . 

75  =  ZE /1000, 

5  II=(I-1)*2+1 

F  2 ( ICl »  II )  =  ?4 

c2  ( IC1  ,  i  I  +1  >  =  75  -  -  -  . - -  •  - 

6  COMTlNUt 

11  PRINT* ,  5  OHT  VP1-  1  TO  Gu  0»  TYPE  0  TO  STOP 
READ*, IFG 
If (IFG, £0.0 ) ST0W 
CALL  I.NiT 
GO  TO  11 
FRO 

SUD  OUT1N.  INIT 
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c 

c 

c 


c 


10 


25 


30 

35 


45 

50 


It  IT  U  lIZ,  T  ION  \OUTI  Ni 

INPUT  JVFkKlUl  LIMITS  AND  EPSILON  (th=*01) 

MAKE  INITIAL  GUESS 

COMMON/ VAI/L  IMi,  Li  M2,  *_'t  ,  IP,  H  ,  N  ,  DEL  ,  Ml  ,Ct . . . 

C0MMGN/MAT/F2C512,4) , A ( b) , A A ( 0 ) , 3 < 51 2 ,9» , f ( c ) , F 3 < 512 , A) 
P‘InT*,30H  TYPE;  IN  LIM1  ,LIM2,  f  I 


A  F AO*»  ulMi,  1.1*2,  EF 
Cf  *1. 

PRINT*  ,  5HLI  H=,L1M1,5HLIM2=,LI  32, 3HEL=,  LI  ,3rtC_=  ,C' 
PRInT*, 25HTYP-  IN  <, G 1 , G2, *  1 , A ’, A  3 
R EAC*, AA 

PRINT*  ,  10  H  1ST  GUESS,  A  A 
C(  L  =  1 

CALL  CORF (At ,  •  J2) 

PRINT*, 7H  .JMlNi=,f  J2 
CALL  L  OOP  ( ,<  J  2  ) 

End 

SUBROUT  tNL  OD 

CALCULATE  C  mAT'-.IX  Fr-.(  M  Q  uNf  R 
D  1  M;  NS  1 0  f  i  R.  (  b  )  ,  Q  (  a  ,  o )  ,  W 1  ( f> ,  o ) 

COMMON /YAk/ LI  -1 ,LIM2 ,Ft ,:°,H,N  ,  JSL ,Nl ,CE 
CDMMC)N/MAT/F2(  512,4)  ,  A  (t)  ,  A  A  ( t  )  ,0(512,9)  ,  C  (  *  )  ,  5  3  (  5 1?  ,  A) 

common/  -tAiNi/x  ):M,uri,-  i 
NOIM=IP_NDI  M1  =  IP  +  1 


DC  10  1=1, IP 

-  ( I )  =  0  . 

:o  io  j=i,i p 

Q(I, JJ =0. 

CO  35  <=N1,N 
SN=F2( <, 3 )  -  F  3  (  K ) 

SM=F2  (<,=*>  -F3 (K, 2) 

CO  30  1  =  1,10 

cc  ?f  j=i ,: p 

Q(I,J)=U(I,J>  *?.*CM*H(K, Z) *3( K 
IF(1.ST. 3.JR.J.3T.3) GC  TO  25 
0 (I, J>  =  0(1, J) *2.*8(<,l+b)*3 (K, 

CONTIN'J: 

R  (I)  =  fs(I)-2.*C  :*SN1*B(K,I) 

IF (I • GT  »  3 ) GO  TO  30 
'  (1)  =  "  (1)  -2 • *  SN*  3 (K  *  1 ♦ 5 ) 

CCNT1NU- 
CCNT1NU- 

CALL  GMi,\V(  IP,  IP,0,vU,MR,l> 
1P1=1°*IP 
DC  50  1  =  1,  IP 
11=1 
0 (I ) =u . 

DC  A5  J  =  I,1P1,IP 
C  (l)=l  (I)  -Ml  (  J  )  *r.  ( II  ) 

11=21*1 
0(I)=-0(I> 

END 

.SUMvOJTINl  LOOP  (  -  J2> 


,J> 
J+t ) 
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C  ITEkATIuN  PROCESS 

C  COMPUTE  A(I  +  1)=A(IU0(I) 

C  DETERMINE  WHEN  MI+1)  IS  acceptable  -  — 

DIMENSION  000(6) 

C  0  MMON  /  V  A  K/  L I  Ml ,  L I  M2 ,  E  E  ,  I P ,  M-,  NvtJtL  *  Nt-fCS - 

C0MM0N/MAT/F2<512,4) ,A<6),AA(6),B(512,9>,D<6> 

NCT =MC T  =  0  - -  -  - 

1  CALL  00 

2  DO  100  1=1, IP  - 

D (I ) =D ( I ) *DEL 

100  A(I>=AA<I>*9(I>  "  - 

IF<A(1>.LT. Q..0R.A<2) .GT.O. )GO  TO  15 

IF(A(4).LT.O..OR.A<5>.LT.Q..OR.A(6>.LTvO.)GG  TO  15  . . 

CALL  COLF (A »RJ1 ) 

6  IF(r  Jl.LT.R  J2)GO  TO  30  — . -  -  -  - 

15  CEL=.5 

NCT=NCT  +  1  - - 

IFINCT.LE.LIMIIGO  TO  2 

PRINT*, 7HLKR0*  1 , 3H  A=  ,  A ,  5H  NCT=,NCT,5H  MCT=,MCT*RJl,fcJ2 - 

PR INT*  *  3H  J=,0, 3HAA= ,«A 

RETURN  •  -  - .  -  - 

30  DO  35  1  =  1, IP 

DOO(I)  = ABS ( 0 ( I ) / A A ( 1 ) ) 

35  IF(0C0(I).GT.lE)G0  TO  20 

GO  TO  40  - - - 

20  MCT=MGT*1 

IF(MCT.LE .LIM2) GO  TO  9  . . * - - 

PRINT*  ,7HER"0R  2,3H  A=,A,5H  NCT=,NCT,5H  MCT= ,MCT ,F J1 ,R J2 

PRINT*, 3H  D=,D,3HAA=,AA -  - 

RETURN 

9  RJ2=FJ1  - — - 

NCT  =  0 

DEL=1.  - - - - — 

DO  25  1  =  1, IP 

25  AA(I)=A(I)  '  - 

GO  TO  1 

40  PRINT 45, (A(I), 1=1, IP) ,NCT,MCT  - 

PRINT*, 3HAA=,A4 

PRINT*,  3H  0=,0  - - - — - 

PRINT* ,5HJMIM=,RJ1 

45  FORMAT (1X,6G12.5,2X,215) - 

END 

SUBROUTINE  COLF(M,RJ)  - - — ~ — 

DIMENSION  M(6> ,IWK( 11) 

COMMON/VAK/LIMl,LIM2,EE,IP,H,N,OCLiNt,C^— - 

COMMON/MAT/ F2( 512,4) , A < 6) , A A ( 6 ) , B ( 51 2 ,9 ) , D ( 6 ) , F 3 C 512 ♦ 4) 

U  =  W(D  - - - - 

W=W(2) 

r=w<3)  - - - - - 

A 1=W ( 4 ) 

A  2=N( 5 )  - - - 

A  3=M( 6) 

_C  =  1 . 3  U  - 
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x  =c*  V U 

X  1=1>C*V 

REWIND  1 

c  compute  partial  derivative  of  mean  error  wrt  xgain 

DC  10  1  =  1, N 

T= ( I-l ) *H  •  -----  - - 

fHI)  =F  2(1,2) 

21=32=100. 

TF(U*T.GE.2Q0.)S1=0. 

IF<C*V*T.LE.-200.)S2=Q. 

IF (Si. NL. 0. >31=EXP(-U*T) 

1F(S2.N,_.0. ) S2=1XP(C*V*T)  -  -  - - 

F  3  ( I)  = (S2*X1-X1-C«C*V« Y*T*S1>  /  (X*C*  V) 

10  F3  ( I)  =  F3  (I)  -(  <U*X1+C*V)*S2-X*XH-C*C*V*Y*S1>/(C'*V*X»X> 

CALL  VCONVOt  F3,t3,N,U,:w<) 

C  COMPUTE  PARTIAL  DERIVATIVE  OF  MEAN  tRROR  WRT  G A  MM A 1 

DO  15  1=1, N 

FF=F3(I)*H  -  - 

WRIT  E ( 1 ) FF 

8 ( I ) =F 2 ( 1 ,2  )  . -  -  " 

T=(I-1)*H 

S 1  =  S  2= 1 0  0  •  -  . 

IF (U*T .GE.2 CO. ) 31=0. 

IF(C"V*T.L£.-200.)S2=0.  - - - 

IF  (SI.  ML.  0.  >  S1=EXP(-U'T> 

IF(S2.Nl.O. )S2  =  tXP(C»V»T)  -  - . . 

F3(I)=C*V*X*(S2«-(  1  +  T*  <  X  +C4U*  Y  )  >-Xl+C*Y*Sl> 

F3(I)  =  F3(II  -<U*2.*C*V>*  (  (U*Xl+C*V)«‘S2-X*Xl*C*0*V*r*Sl) 
15  F3(I)=F3(1) /(C*V* V»X*X) 

CALL  VC ON VO (F3»3»N,N,lWtO  - - 

C  COMMUTE  PARTIAL  DERIVATIVE  OF  MEAN  E  RRO^  Wr.T  GAMMA1 

DC  17  1=1, N 
FF=F3(I»*H 

WFITEdIFF  . . - 

B(I)=F2(I,2) 

T=(I-1)*H 

S1=S2=100. 

IF(U*T.Gf .200. ) S1=0.  -  - - - 

IP(C*V*T. LE.-200. )S2=0. 

IFCSl.ML. 0. » S1=EXP(-U*T>  "  -  " 

IF(S2.Mfc.0.)S2=£XP(C*V*T) 

17  F3(I)  =  (U*S2 -U-C* V*C*V*  S 1 > / ( V*  XT - 

CALL  VCUNVO ( F3 , 6 , N, N , I WK) 

COMPUTE  ESTIMATE  OF  OT  —  - . . . ~ . 

00  20  1  =  1, N 

FF=F3(I)*H  -  — 

WRITE( 1) FF 

B (I>  =F2 ( 1 ,2  )  -  - 

T= (1-1 ) *h 
S1=0. 

IF(UM  .LT.2  00.  >S1  =  EXP(-U*T> 

20  F3 ( I ) =S 1 

JIALL  VCONVO(F3,R,N,N,IWK) 
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8  1 1 1 3 ) =F2(1,2> 


1C  COMPUTE.  PAt.TI>.L  OF  OT  ESTIMATE  HPT  KGAIN 

00  25  1  =  1, N 
T  =  <1-1 >  *H 
FF=-F3(I)*H*F2(I) 

:  8(1,7) =  l  • 

!  B (I, 8) =FF  . . . .  ■ 

S1  =  0. 

1F(U*T.LT.20Q. >$1=EXP(-U*T) 

F  3( I) =T*S1 
B  ( I) =F2(I ,2) 

IF(1.E0.1)G0  TO  25 

3  (I,*)=  (Bd  ,8) -8(1-1  ,b)  )/H  .  . .  . 

25  CONTINUE 

CALL  VCONVO (F3,3,N,N,IWK) 

B (1,6) =0. 

!  no  27  1=1, N  -  -  - 

3(1,5)  =F3(I)+H 

IFd.EQ.DGO  TO  27  - - - - - 

B(I,P)=(B(I,5)-B(1-1,5))/H 
27  CONT I N  UE 

C  COMPUTE  PARTIAL  DERIVATIVES  OF  VA-IANCE  WET  ALPHA  ( J) 

DO  40  J= 1 , 3 
DO  35  1=1, N 

T=  ( 1-1 )  *  H  . .  . .  . . . . — ' 

II=J*b 

3(I)=‘3(l,lI)#b(I,II) 

S 1  =  S2= 100  • 

IF(U*T.GT.2n0)Sl=0. 

IF(C*V*T.LT  ,-2QU.)S2=0. 

lF(Sl.Nt  ,  0.  >  S1  =  E  XP(-U*T>  -  -  -  - 

IF(32.NE. 0.) S2=EXP(C*V«T) 

F3C  >  =  (C*  <S2MX+C*U*Y) -U*C*Y*S1)/X) ♦♦2 
35  CONTINUE 

CALL  VCONVO (F3, J,N,N,IWK> 

DO  37  1=1, N 

FF  =  F3(I)*H  "  •"*  - - * 

37  WAITE(1)FF 

40  CCNTINU.. 

C  COMPUTE  PARTIAL  DERIVATIVE  OF  VARIANCE  WRT  KGI A  N 

DO  50  1=1, N 
T=(I-1)*H 

S1=C2= 100 ,  -  - 

1F(U*T,GT ,200)31=0, 

IF(C*"V*T  ,LT  ,-200  .)S2=0, 

Ic(Sl.Ni, 0. ) Sl=uXP(-U*T) 

|  1  lF(S2,Nt.0. >S2=EXP(C*V*T> 

|  5  (1,7) =S2*(U*C*Y+X)-U*C*Y*S1 

I  B  <I,4) =Al*A2*B(:,c)*B(I,8)»43*3(I,9)*3(I,S) 

f  6d)=0(i,4) 

50  FJ(l)  =  2.*C*C*3<:,7)MXMr2*Xl-C*r*SlMl*tJ*T) ) -312 ,7T1  r(X*r+XT  — 
CALL  VCONVO(F3,3,f;,N,IW<) 

.4.  w  _DO  c3  I  =  1 , N 

y 


? 
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53 

C 


55 

C 


60 

C 


65 


75 

77 


T  =  ( I  - 1 )  ♦  H 

9(1*3) =F 3(1) *H  . - . 

6(I)=A2*B(I*'>)*£(I*5)+A3*8(I,'j)*3(I,6) 

F  3 (I ) =2«*C*C*B(I ,  7)  ♦  !3  (I *7)/ (X*X)  - - -  - 

CALL  VCONVO< F3,J,N,N,:WK) 

coipur:-:  partial  derivative  of  variance- wf-t  (tamm-ai - 

DC  55  I=1,N 

FF=F3d)*H*8(I,3)  - - 

Wf.ITE<  1)  FF 

T  =  ( I -1  )  ♦  H  -  -  - 

Sl=100 

IF(C*V*T.LT.-20Q.)S1=0.  . . .  " . 

IF(Sl.Nt.  0.  )  S1  =  EXF(C*V*T) 

9  (I)  =  Ml,  L)  '  -  - - 

F3  (1)  =C*C*C  *2,  *3  (1 ,7  )M  Si*  <  T*  (  X«-tJ*C*  Y  )  +1 )  *  X-0  ( I  ,7>>/(X*X*X) 

continue  .  - . 

CALL  VCONVD(  F3,  i.N.N.IWO 

CCMPjr'T  PhKTI^L  u!  RIVF  TTV*  OF  VARIANCE  MPT  'GAMMA  2 
DC  60  1  =  1,  N 

T  =  (1-1 )  *h  - -  — . .  . 

FF=F3( I >*H 
WRIT-  (1IFF 
St  =  S2= 10 0  . 

IF (U*T  «oT .200. I jl  =  0. 

IF(C«V*T.lT  .-210. >S2=C. 

I F (31 •  NL.  0.  )  SI* X°( *U*T  )  -  ”  -  — -  . 

lF(S2iMF.  0.  )  S2  =  :  XP(C*  V*T  ) 

B(I)  =  3  (!,*•) 

F3(I)=C*C*C*2.*U#B(I  * '/  )*(32-Sl>/(v*X) 

CONTI  NU :L 

CALL  V  C  ON  VO  (  F  3  *  -j  ,  N  ,  N  *  I W  <  ) 

COMMUTE  MlaN  T- ACTING  'KJ1  . . .  '  “  . . .  . 

10  6‘.:  I=1,N 

T=(I-1)*H 
FF  =  F3 (I  )*H 
WFIT-.(  1 )  FF 
R(I)=F 2(1*2) 

S1=32=1PC.  - . .  '  . 


IFOJ+T  ,&T  .200.  )  i"l  =  0. 

IF (C*V*T • LT .-200. >$2=0. 

if (si. Mr • o. )Si=txp<-u*T) 

IF(S2.Ni,0.  > 32=iXF<0*V*T> 

F3(I)  =  (  (U*Xl*C*V)*S2-X*xi+C*C'1  V*Y'*S1)/(C*V*X>  f* 

CONTINUE 

CALL  VCONVO( F3.3»N»K*iW<) 

REWIND  1 


DO  75  1=7,1 
DC  75  J  =  1  *  N 
R  E  AC  (1H(J,II 
DC  77  I  =  a  ,6 
DO  77  J=1,N 
i AO  ( 1 )  i( J,I  ) 
DC  .0  1=1.3 


piilt  IS 


3tfiT  «  - 


********** 

,  .nfl 


72 


f\J  fO 

«o  «o  o 


so 


’W 


3 


c 


85 


CO  .a0  J=1,N 
READ (1)B(J«I) 

RJ=0  * - - 

COMPUTE  STANOARG  DEVIATION  OF  TRACKING  ERROR 

00  83  1  =  1, N  "  '  -  - 

F3II,2)  =SQRT(3(I,4)*AH-B(I,5>*  A2  +  BC1 ,6)*A3) 

DO  82  J=l,b  - — 

B(I,J>=.5*3 (I, J) /F3< 1,2) 

CONTINUE 

CALUCATt  VALUE  OF  COST  FUNCTION 

DO  85  I  =  N  1 ,  N  “ . .  ~ . — 

F  3  ( I)  =F3  ( I)  *H 

RJ=RJ*(F2  Cl  ,3)-F3(I>  )  **  2  +  CE*  C  F 2  C I ,  L)  -T3tT,^TT**-2  ~ 
END 

SUBROUTINE  GMINV  <NP.,NC,  A  ,U,‘MP  ,MT) - 

DIMENSION  A<1),U(1),S(10),IP(10) 

COMMON/ MAIN  1/  NOIM.NOIMl  -  — - -  - 

T0L=l.t-12 

ADV=T0L*T0L  . . . . . 

MP=NC 

NRMisNft-1  - 

NCMl  =  N  C-l 
JJ=1 

00  5  J  =  1 , NC 

S(  J)  =  D0T (NR, AC JJ) ,A( JJ>)  -  - 

5  JJ= JJ+NDIM 

11  =  1  - - 

DO  19  1=1 ,NCMl 

XMAX  =  SCI)  -  - - 

I MAX=I 

IFCNC.EQ.l)  GO  TO  15  - 

IP1=I+1 

DO  10  J  =  I  PI ,  NC  - - 

IFCXMAX.GE.S ( J> )  GO  TC  10 

IMAX  =  J  -  - - 

X MAX  =  S  C  J) 

10  CONTINUE  "  - 

IFCIMAX.EQ.I)  GO  TO  15 

SCIMAXI=SCI»  - - 

SCI) =XMAX 

KK=(IMAX-1> *NDIM*1  •  - -  - - - 

CALL  SNAP  (ACKK)  « A ( 1 1 ) , NR, 1 ) 

15  IP(I)  =  IMAX  - - -  “ 

19  II=II+NGIM 

T  0L1=S ( 1 ) ♦AOV  -  - 

A0V=T0L1 

JJ=1  .  . •' - 

DO  100  J=l, NC 

F  AC=S  C  J)  - 

JM1= J- 1 

JPM= JJ ♦NPMl  - - 

JCMaJJ+JMl 
DO  20  I*JJ, JCM 
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-•  -i*t  £A  U*  •  ** J Mm***-* 


20  UCI>=0. 

U  (JCM) =1.0 
1 F ( J. EQ. 1 )  GO  TO  5<* 

KK=1 

00  30  K=1,JM1 
I F (3 ( < ) • EQ. 1.0)  GO  TO  30 
TEMP=-DOT  CNR,A<  JJ)  ,A  (KtO  > 

CALL  VAOO(K,TLMP,U(JJ> ,U(K<) ) 

30  KK=KK*NDIM 
00  50  L=l,2 
KK=1 

DO  50  K=l,JMi 
IF(S(K)  .EQ.  0.  >  GO  TO  50 
TEMP=- JOT  (NiR,A  <  JJ)  ,A  (KtO  > 

CALL  tfACO(N= ,TEMP,A( JJ) ,A(KK) ) 

CALL  VADD(K»TEMP,U(JJ)  ,  U  ( K< )  ) 

50  KK=KK+NOIM 

TOLl=TOL*FAC«-AO\f 
FAC  =  COT  CNF.*  A  (JJ)  ,A(JJ>  ) 

54  IF(FAC.GT.TOLl)  GO  TO  70 
DO  55  I=JJ, JRM 

55  A ( I ) =  0  • 

S ( J)=0. 

KK=1 

DO  65  K=1,JM1 
I F (S ( K ) ■ EQ. 0  ■  )  GO  TO  65 
TEMP=-DOT(<  ,'J(KK)  ,U(  JJ)  ) 

CALL  V ADO (NR, TEMP, A  (  JJ)  ,A (KK) T - 

65  KK=KK+NOIM 

FAC=D3  T ( J,U ( JJ) ,U( JJ) ) 

MR=M*-1 
GC  TO  7 5 
70  S(J)=1.0 

KK=1  . -  - 

DC  72  K=l» JM1 
I F (S ( K ) .FQ. 1. )  GO  TO  72 
TEMP=-DOT(NR, A( JJ) ,A (KK)  > 

CALL  VADO(K,T£MP,U(JJ) ,U(KK) ) 

72  KK=KK+NDIM 

75  FAC=1. /SORT (FAC)  - 

DO  80  I=JJ,JPM 
80  A(I)=A(I)*FAC 
DO  85  1  =  J J  ,  JCM 
85  ll(I)*UIII*FAC 
100  J J=  JJ+NDIM 

I F  ( MR .  EQ . NP  .  OF  .  MR .  ET) . NC ) '“GO  “TCT ITT 
IF(MT.ME.O) PRINT110,NR,NC,MP 

110  F  ORMAT  (I3»1HX»I2»8H  Mt  RANK,I2Q - 

120  NEND=NC*NDIM 

J  J=  1  . . 

00  135  J=l»  NC 

00  125  1*1,  MR  -  .  . . 

1I=I-J* JJ 
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125  S  <I»  =  30T2  <N-NO,U(JJ>  ,A  (II>) 

N  LN j=N  EN  L“N  0 1 M 
11  =  J 

00  130  1=1, NR 

J (II ) =  3 ( I )  -  - 

130  II=1I+NLIM 
135  J J= JJ+NOI Ml 
00  33  J=2,NG 
K=NC-J+1 
KK=IP( K) 

1F(K.NE.1P(K>)  GALL  SWAP  (U(  K  )  ,  'J(  KKIvNT,  NCIM1 
39  CONTINUE 
S  ETU'.N 
PNr 

SUBROUTINE  VADOCNfCl ,/ ,3) 

DIMENSION  A (1) ,0(1) 

00  1  1=1, N 
1  «(I)=A(I)+C 1*3(1) 

RETURN 
£  NO 

^UNCTION  CO  T  ( N  E ,  A  »  H ) 

DOUBLE  PRECISION  DUTIABLE 
0 I Mb NS I ON  A ( 1 ) ,7(1) 

OOT1=0.do 

IF(riK.LE.O)  GO  TO  2 
DO  1  1=1, NR 

1  00Tl  =  0UTH-0BLi.  (A  (I>*  B(I>  > 

2  OCT=DOT  1 
RETURN 
END 

SUBROUTINE  SW*.P  (  A  ,B»U»  INC) 

DIMtNSION  A ( 1 ) ,3(1) 

NN=N*INb 
1  =  1 

1  IF(I,G7«NN)  RETURN 
T  F MP  =  A ( I ) 

A ( I) =6 ( I) 

3 (I)=TEMP 
1=I+1NC 
GO  TO  1 
LN'O 

FUNGI  ION  00  T2 ( NN , A , 3 ) 

DCUBLt  PRECISION  t'OT  2  ,C'6LE 
DIMENSION  A(i) ,3(1) 

COMMQN/MA INI /NOIM 
OOT2=0 ,00 

IF (NN,  Lt»  0)  GO  TO  2 
DO  1  1=1, NN, NOIM 

1  C'JT2=D01  2*03LE(  A  <I)*B(I>  ) 

2  D0T2=00T2 

TU-N 
F  NO 
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2 

1 


100 

200 

.01 

5. 

-2. 

*-4 

• 

• 

• 

CNJ 

1 

•  1 

1 

100 

20  C 

.01 

10. 

-6. 

-6.  .001 

.01  .01 

1 

10  0 

200 

•  01 

2. 

-10, 

-D.  .001 

.01  .01 

1 

100 

200 

.01 

1* 

-1. 

-1.  .001  . 

001  .001 

1 

100 

200 

.01 

1. 

—4  • 

—4.  .0001 

.0001  .0001 

0 
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APPENDIX  D 


PROGRAM  LISTING  OF  COMPUTER  SIMULATION  OF  AAA  TRACKING  TASK 


COMMON/  VAF/DEL»N»M»CT0D»ECC»ET,T»C7O 

'PPtN'T*,  20HUSE  eveKy - SAmpLCST 

REMIND  Z 


REAO*,Il 
C70*fl. 

CEL*I1*»0 23 

T5iwn - 

N*3$N1*9 

*>EAO*,IFG 

T2*i 

IF  (IFG. EC. 2)12=2 


TS7,  “^FCfTET  6ft  rrCR  “EOtF 


Pmsl.Z8IP<y,2>  =  1.3<* 

P  <  8 J  *p ( !,2)=0. 

DO  15  IC=l2f2 

PRIhT*,«iOpTYFE  IN  GAMMA1»GANK  \2,KGAIN,Pi,F2,P3  FCf-,lC 

IT  TfEJnF*<MJ,IET73*T»T? 

DO  105  1*1,1050 

- READ  12,314,3,  CTcTd'TeTVC'.eCC - 

_  KK*I»I1«1  _ _ 

ICi*ICl4l* 

- t=-tTC  1  -IT  *DFI . . . . 

CALL  CBSEL(ELERR,S0EL,C7,P(9>) 

- IF  (IFG.NE.2) CALL  CeSAZ<AZERP,SCAZ,C7,P> - 

WPITE  (6  ,* )T , AZERR , ELEPR  «SQAZ , SSEL 
I  F  (MCDTI ,  : 


1W 

3 

5 - 


> ,  t  »  n  tc''',tkcnr(»?TjiZ  ,SDEL 
IFCIC1.EC,1)FRINT5,T,AZERR,ELERR,SDAZ,SDEL 
CONTINUE 
PORMAT (6C12.5) 

FORMAT (5G10*  3J - 

END 

SUBROUTINE  OESEL  rELERP»5CEL  »C7~,T1 - 

COMMON/ YAF/DEL,N,Nl,OTOD,EOO,ETtT 

DIMENSION  A(3,3),Z(3)*XX(3, 3)",Rl  (3)  ,Kin,3»  ,N2<3,3>  ,TTf> 
OIMENSICN  F2<3> ,GAM1P(3) 

IF (T* GT..501 ) GO  TC~7 - 

zm*-P<7) 

^ZTZJTZTJnvTJT *rm - - 


DO  10 
HTJT- 


j*l ,N1 

IT, - 


10 


XX(J)=0. 

ni)«P(3)4P(7)»  (P(i)»pm»p<gn - 

A<2»*-PC3l*PI3)-PI7»*PC3)*fPCl)*PC3)*P12n 

■mTiimTFFTZI - - 

A<fU-Pm-P<7l*P«3l*P<2) 

H7)«-PI7»*PI2I - 
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“J<8)*P(7>«P(3)»P<2> 

- riwn] - 

call  DSCR1(N,A,CEL,W1,W2,1C) 

F2Tir=“PT7) 

F2 (2>=-PC7)+P(3> 

- r2TJ)=0% 

00  36  K=1,N 

Trmrm~zr. - 

00  36  J  =  1,N 

~~  36  (TANiF  (TO^MIPTKJ +W?TK7J1*FZ(  J  > 

TO  50  1  =  1, N 

ii  *1  . ‘ 

Pi(I>=0. 

00  45  J3I,Ni,N 

_  pici)=p.i<n»w2(j>»z  (ii) _ _ 

45  11*11+1 

50  CONTINUE 

00  60  X  3 1 ,  N 

_  7(I)=«?» _ 

00  60  JslyN 

60  A <I,j»=*l(I>*Ri (j) 

00  20~I3i,N 

__  _  I3=l+N _ 

I2*I3+N 

2_0 _ Pl(l)*M2(l3) +H2(I2) _ 

7(l)3Z(2)sZ(3)aO. 

xxd)=.oeci _ _ _ 

1  P4aZ(2)-Z(3l*P(3:>'Z(i) 

P5=(P4-P(6))/C£t 
IF  (T«  EQ«0 • ) PSscCO 

V= (P ( 4) +P  (5)«P4»P4+P(6)*P5»PS)/CEL _ 

P(8)*P4 
CO  25  1=1, N 
II  si 

M2(I)=0. _ _ 

DO  15  J*I,N1»N 

wzdiswzg)  ♦wi(j)*zcin _ 

15  11*11+1 

25 _ continue _ 

DC  35  i*l,N 
Z(I)=N2(I)+R1(I)*EDQ 
35  CONTINUE 

UBAR=-P (1)*Z(1)»P(2)*(Z(2)+P(3)*Z(1)«Z(3)) 
IF (ABS(U6AR) .LE..785)  GO  TO  3t 
UFIXs.7€5 

tF(L'BAR.<ST.O.)  UFIXS-.765 
DO  37  J=l,N 

37  Z(J)*Z(J)+GAH1P  <J)+ (UBAR+UFIX) 

30  CONTINUE 

CALL  MULT(Nl,XX,N,NlyH2,10) 

00  40  1*1, Ni 
Of  XX  (I) *A  IX  )*V+H2  II > 

__30EL*SQPt  (XX  (1)  )  +  150/3. 14156 
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C7=£T-Z<1> 

EL ERR=Z(i)*18P./3. 14156 

End 

SUBROUTINE  CESA  Z  <AZEF!i,S0AZ,C7  ,P> 

Tl  KENSlCN  "A A  ( 3 , 3  f,  ZI  <  3  >  ,  *<3, 7 >  ,p2~ ( 3  )  ,  WK  ( 3  >  ,  W  3  ( 3  ,  3  >  ,  k A  (  3',  3 ) 

_ OIMENSICN  GAH1P(3) ,P(6) ,F2<3) 

COMMON/ VAF/O EL,  k,M,CTCr ♦  EDO , ET , T ,C7C 

IF  CT« 6T«« Ofll )GO  7 C  1 _ _ 

R2ti>=-F<7> 

R2(2)=R2(3)=F  <7>*P(3> 

DO  10  JsltNl 
AA ( j)=0. 

10  ><J)=0. 

Z1C1)  =  Z1<2)  =  Z1(3)=?« _ 

X(i)=. 01005  *  '  . 

Zl (3>=0. 

1 - F<*=~ZTT2  >-  Zl UT*WU )  *  ZI  111  ~ 

P5=(P4-P(8))/DEL 

IF(T.EG.3.)P5=07CC 

V= (p ( 4 ) ♦P(5)*PA*PA+P(6)*F5*°5)/tEL 

P<8)=P4 

070= ( C7-C70) /DEL 
ETiiC0S<C7) 

02  =-T  AN(C7)*C7C 
C7C=C7 

AA  <l)=02*Cl*P<3Mni«pm*<P<l)*P(2)*F<3> ) 

AA«2)=-D1*P(3)*P(3)*P(3)*C2-P(7)*F(3**D1*CP(1)*P(2)*F(3) ) 
AA(4)=D1*C1*F(7)*F  2) 

- - AA i 5) = -C I*¥  f  2) - F (7  »^P (2 ) * P  C 3 ) «D1 

AA C7)=-P<7)*Ll*P(2) 

AA  <8>  «STT7T*P  (Z)  *F  T37  *01 
AA (s)=-P(3)*Ci 

CALL  0SCPY(N,AA,DEL.N3,W4»ir) 

00  25  1=1, N 
Il^i 

w< (i) =o« 

00  15  J=I »N1 , N 
HK(X)*MK(I>*M3(J)*Z1(XX> 
lyilTfl 
continue 

n0-J5~T*T,H  ~  ‘ 

71(I)=H«(I)MH<i  CX  +  3)  ♦WMI+61  )  *0l  DD 
CONTINUE  ' 

UBAR=-P<1)*Z1  Cl)-F(2>*CZU2)*PC3)*Zl(l)-ZiC3)l 
IF(ABS<ueiR)  .LE.,  829)  GO  TO  78 
F2 (1)=P  171*01 
F2C2)«“F(7)*P(3)*C1 
F2< 3)=0. 

00  3f  k*I7N  " 


6AH1P(N»*8. 


36"^AH1P(<)=CAH1P(k) *W4 ( < , j ) *F2 <j ) 

ufixs.bm 

TF<UBAR.GT.O.)  UFlX=-.829 
CO  37  JsltN 

37  7i  C  j)=Zl(J>  ♦CA*UP<j)*<UBAS  +  UFIX) 

38  COKflNUE 

HO  50  X=1,N 

II  *1  '  . . . . 

WK ( I) =0 

DO  4ff  J=I  *  Ni  » N 
WK(II=WK(I)  +W4( J)  *‘-{2  (III  «Di 
40  XI»XX*1 

50  CONTI MU E 

-  “CALI  HUlTiW3,X7F,Ni,*4,iO> 

DO  6(5  I  =  1,N 

do"6(T“j=L7^  ~~  ~ . . 

II=<j-i)«K+I 

60  jrmrsWKfrnNKt  j)*v*w4(H> 

?OAZ=SQRT <X<1> )  *18!)./<3.14156*C1> 
AIEMsT" nU*  1 8T  .  /T3.T  4 15  fi  i T 

END 

- su  eko  uti  Nr  ^urnE  ,  f ,  l  ru;*,  m  ~ 

DIMENSION  E  (LI)  ,F  (LU  ,G(9)  ,H(H) 

- D'OT'i'f  'T=iTl 

11=1 

OCT1C  T=IYL 
tEMF=T. 

00  5  j=i,Ti»L 
tEMP=TENP*E(j)*F<II) 

5 - :?=Tl*i - 

JCK-JK-1 >*L+I _ 

”  H  IKK)  =TeK? 

10  G<<K)=TEMF 

IF  (NR.EC,TTffETUT?N 
DO  20  I=1»L 

co  20  K=IfL 
temf=o. 

. . .  TI  =  K . .  . 

PO  15  JsX'Ll'L 
“  fEMP=TEFF4GCJ)*mi) 

15  TI  =  I1*L 

KK  =  (K-l  )*L  +  I 

23  HI-fK^TCNP  _  _ _ 

12=1-1 

CO  30  1=1,12 _  _ _ 

L3*I*1 

_ ro  30  J  =L 3 1 L _ 

K1«CI-1)*L4J 

_ F2L=<J-l>»l»I _ 

30  H(K1)=H(K2) 


END 
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Tu  8  R  C  U  T 1 N  E  b <5 c R T  { K 0 1 M ,  A  , C E L  ♦  E  A  *  E A I N >  , NT > 

DIMENSION  A  tl),EA  <1 > «EAlNT < i> . CC£F < 3 0_> _ 

SttS  EA»E*P(A*CEL>,eAINT=TNTEGRAL  EA  C  TC  DEL 

NCIM1=NCIK*1 _ _  _  _  ___ _ 

NNsKDIH^KDIM 

NTM=NT-1  _  _ _ _  _  _  _ 

CCEF(NT>=1. 

PO  10  I=1»NTH1 _ 

lisNT-I 

PQEFC II  >  =CEL*cCEF  <!!♦!  I/F LO  AT  { l) _ _ 

NT  MUST  BE  AT  LEAST  3 

.CALL.  DIAGCNJIK.EAINT.A.CCEFC1).C0EF(2))  _____ 
TO  60  L*3,NT 

call  HULTCAtEAINT ,ncim,nn,ea,i )_ _ 

IF  (L.EC.NT)GC  TO  70 

CALL  PIAG  (NOlNt  EA 2l K T, E  A  jjL  •  Pi  C OEFJi  LL_  _ 

DO  90  Ilsi,NK,NCm 

EA  ( II )  =  EA  (II  >  »1  «C _ _ _ _  _  _ 

CONTINUE 

END _ _ 

SUBROUTINE  DIAG <NDIM , A  ,B ,C1 *6 2 ) 

DIMENSION  A  (1 )  t  E  (1 )  _____  _  __  _ 

NDIN1=NCIM+1 

_  n_n=ndj.n*ncim _ _  __  _  _  __ 

NM1=NDT>-1 

ti»i _ _ 

IF (Cl • E  T. 1 . 0 )  GC  T 0  1C 

PO  5  J=1  «NN» NDIK _ _ _ _  _ 

K=J+NM1 

CO _ 4  I»J»K _ _  _  _  __ 

►  A (I )=Ci *E (I) 

A  (II)  -  A  <  1 1 )  +C2 _ _ 

;  u=Ti+NcI>T 

»ETUf<N  _  _ _ 

i  ro  7  j®l,KN,NLIN 

_  K^J*N_Mi _ _ _ _  _ _ 

DO  6  I®J»K 

>  A(I)®B(I) _ 

A ( II) *A (II)  *C2 

'  n=n*Ncm _ _ _  _  _  _ _ 

RETURN 

END 
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